#pragma once

void compute_OPC(O5_ccccssss OpmPC, gamma_ccss GOS, gamma_ccss GTW, int ipm);
void compute_OPV(O5_ccccssss OpmPV, gamma_ccss GOS, gamma_ccss GTW, int ipm);
void compute_OpPC(O5_ccccssss OpmPC_out, gamma_ccss GOS, gamma_ccss GTW);
void compute_OmPC(O5_ccccssss OpmPC_out, gamma_ccss GOS, gamma_ccss GTW);
void compute_OpPV(O5_ccccssss OpmPV_out, gamma_ccss GOS, gamma_ccss GTW);
void compute_OmPV(O5_ccccssss OpmPV_out, gamma_ccss GOS, gamma_ccss GTW);
void compute_OPC_2r(r_O5_ccccssss OpmPC_out, gamma_rr_ccss G,int ipm);
void compute_OPV_2r(r_O5_ccccssss OpmPV_out, gamma_rr_ccss G, int ipm);
void compute_OpPC_2r(r_O5_ccccssss OpmPC_out, gamma_rr_ccss G);
void compute_OmPC_2r(r_O5_ccccssss OpmPC_out, gamma_rr_ccss G);
void compute_OpPV_2r(r_O5_ccccssss OpmPV_out, gamma_rr_ccss G);
void compute_OmPV_2r(r_O5_ccccssss OpmPV_out, gamma_rr_ccss G);


/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// O_g1g2^{+-} (0)= 1/2 [ bar{psi1}(0) g1 psi2(0) bar{psi3}(0) g2 psi4 +-  bar{psi1}(0) g1 psi4(0) bar{psi3}(0) g2 psi2 ]  (eq 26 bibbia)
// -> G_g1g2^{+-}= 1/2 [ GOS*GTW -+ GTW*GOS ]
// note the change in the sign +- due to the fact that psi are fermionic fields and thus grassmann variables
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

///////////////////////////////////////////////
/// Note that:
/// G  ~  GOS*GTW +- GTW*GOS 
/// with GOS=G_00 and GTW=G_01 or GOS=G_11 and GTW=G_10
/// is equivalent to
/// G ~ G_r1r2*G_r3r4 +- G_r1r4*G_r3r2
/// if r=00,01 or r=11,10
/// but not in other cases. For example r=01,11 -> G~G_r1r2*G_r3r4 +- G_r1r4*G_r3r2 ~ G_01*G_11 +- G_01*G_11 ~ GTW*GOS +- GTW*GOS
/// However the only relevant combinations are 0001 & 1110 (the others are equal) so one can use the following formulas
/// For other r-combinations (1000 0100 0010  1101 1011 0111) we can use the formulas in *_2r routines
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////


void compute_OPC(O5_ccccssss OpmPC, gamma_ccss GOS, gamma_ccss GTW, int ipm)
{

	
     // 1. Initialize
      for (int ic1=0; ic1<3 ; ic1++) for (int ic2=0; ic2<3 ; ic2++) for (int ic3=0; ic3<3 ; ic3++) for (int ic4=0; ic4<3 ; ic4++)
      for (int is1=0; is1<4 ; is1++) for (int is2=0; is2<4 ; is2++) for (int is3=0; is3<4 ; is3++) for (int is4=0; is4<4 ; is4++)
      for (int iop=0; iop<5; iop++)
      {
                OpmPC[iop][ic1][ic2][ic3][ic4][is1][is2][is3][is4]=0;
      }

        // 2. Compute for each r combination with one r different


        for (int ic1=0; ic1<3 ; ic1++) for (int ic2=0; ic2<3 ; ic2++) for (int ic3=0; ic3<3 ; ic3++) for (int ic4=0; ic4<3 ; ic4++)
        for (int is1=0; is1<4 ; is1++) for (int is2=0; is2<4 ; is2++) for (int is3=0; is3<4 ; is3++) for (int is4=0; is4<4 ; is4++)
        {


         // VV+AA

         for (int imu=1; imu<5 ; imu++)
         {
                int imu5=imu+11;
                OpmPC[0][ic1][ic2][ic3][ic4][is1][is2][is3][is4]=OpmPC[0][ic1][ic2][ic3][ic4][is1][is2][is3][is4]+
                1./2.*(GOS[imu][ic1][ic2][is1][is2]*GTW[imu][ic3][ic4][is3][is4]+GOS[imu5][ic1][ic2][is1][is2]*GTW[imu5][ic3][ic4][is3][is4])
                -pm(ipm)*1./2.*(GTW[imu][ic1][ic4][is1][is4]*GOS[imu][ic3][ic2][is3][is2]+GTW[imu5][ic1][ic4][is1][is4]*GOS[imu5][ic3][ic2][is3][is2]);
        }
         // VV-AA

         for (int imu=1; imu<5 ; imu++)
         {
                int imu5=imu+11;
                OpmPC[1][ic1][ic2][ic3][ic4][is1][is2][is3][is4]=OpmPC[1][ic1][ic2][ic3][ic4][is1][is2][is3][is4]+
                1./2.*(GOS[imu][ic1][ic2][is1][is2]*GTW[imu][ic3][ic4][is3][is4]-GOS[imu5][ic1][ic2][is1][is2]*GTW[imu5][ic3][ic4][is3][is4])
                -pm(ipm)*1./2.*(GTW[imu][ic1][ic4][is1][is4]*GOS[imu][ic3][ic2][is3][is2]-GTW[imu5][ic1][ic4][is1][is4]*GOS[imu5][ic3][ic2][is3][is2]);
        }

        // SS - PP
                OpmPC[2][ic1][ic2][ic3][ic4][is1][is2][is3][is4]=OpmPC[2][ic1][ic2][ic3][ic4][is1][is2][is3][is4]+
                1./2.*(GOS[0][ic1][ic2][is1][is2]*GTW[0][ic3][ic4][is3][is4]-GOS[5][ic1][ic2][is1][is2]*GTW[5][ic3][ic4][is3][is4])
                -pm(ipm)*1./2.*(GTW[0][ic1][ic4][is1][is4]*GOS[0][ic3][ic2][is3][is2]-GTW[5][ic1][ic4][is1][is4]*GOS[5][ic3][ic2][is3][is2]);


         // SS + PP
                OpmPC[3][ic1][ic2][ic3][ic4][is1][is2][is3][is4]=OpmPC[3][ic1][ic2][ic3][ic4][is1][is2][is3][is4]+
                1./2.*(GOS[0][ic1][ic2][is1][is2]*GTW[0][ic3][ic4][is3][is4]+GOS[5][ic1][ic2][is1][is2]*GTW[5][ic3][ic4][is3][is4])
                -pm(ipm)*1./2.*(GTW[0][ic1][ic4][is1][is4]*GOS[0][ic3][ic2][is3][is2]+GTW[5][ic1][ic4][is1][is4]*GOS[5][ic3][ic2][is3][is2]);

        // TT
        for (int isigma=6; isigma<12; isigma++)
        {
                OpmPC[4][ic1][ic2][ic3][ic4][is1][is2][is3][is4]=OpmPC[4][ic1][ic2][ic3][ic4][is1][is2][is3][is4]+
                1./2.*(GOS[isigma][ic1][ic2][is1][is2]*GTW[isigma][ic3][ic4][is3][is4])
                -pm(ipm)*1./2.*(GTW[isigma][ic1][ic4][is1][is4]*GOS[isigma][ic3][ic2][is3][is2]);

        }


        }// end ccccssss 


}
void compute_OPV(O5_ccccssss OpmPV, gamma_ccss GOS, gamma_ccss GTW, int ipm)

{


     // 1. Initialize
      for (int ic1=0; ic1<3 ; ic1++) for (int ic2=0; ic2<3 ; ic2++) for (int ic3=0; ic3<3 ; ic3++) for (int ic4=0; ic4<3 ; ic4++)
      for (int is1=0; is1<4 ; is1++) for (int is2=0; is2<4 ; is2++) for (int is3=0; is3<4 ; is3++) for (int is4=0; is4<4 ; is4++)
      for (int iop=0; iop<5; iop++)
      {
             OpmPV[iop][ic1][ic2][ic3][ic4][is1][is2][is3][is4]=0;
      }



        for (int ic1=0; ic1<3 ; ic1++) for (int ic2=0; ic2<3 ; ic2++) for (int ic3=0; ic3<3 ; ic3++) for (int ic4=0; ic4<3 ; ic4++)
        for (int is1=0; is1<4 ; is1++) for (int is2=0; is2<4 ; is2++) for (int is3=0; is3<4 ; is3++) for (int is4=0; is4<4 ; is4++)

        {
        // VA+AV
         for (int imu=1; imu<5 ; imu++)
         {
                int imu5=imu+11;
                OpmPV[0][ic1][ic2][ic3][ic4][is1][is2][is3][is4]=OpmPV[0][ic1][ic2][ic3][ic4][is1][is2][is3][is4]+
                1./2.*(GOS[imu][ic1][ic2][is1][is2]*GTW[imu5][ic3][ic4][is3][is4]+GOS[imu5][ic1][ic2][is1][is2]*GTW[imu][ic3][ic4][is3][is4])
                -pm(ipm)*1./2.*(GTW[imu][ic1][ic4][is1][is4]*GOS[imu5][ic3][ic2][is3][is2]+GTW[imu5][ic1][ic4][is1][is4]*GOS[imu][ic3][ic2][is3][is2]);
        }

        // VA-AV
         for (int imu=1; imu<5 ; imu++)
         {
                int imu5=imu+11;
                OpmPV[1][ic1][ic2][ic3][ic4][is1][is2][is3][is4]=OpmPV[1][ic1][ic2][ic3][ic4][is1][is2][is3][is4]+
                1./2.*(GOS[imu][ic1][ic2][is1][is2]*GTW[imu5][ic3][ic4][is3][is4]-GOS[imu5][ic1][ic2][is1][is2]*GTW[imu][ic3][ic4][is3][is4])
                -pm(ipm)*1./2.*(GTW[imu][ic1][ic4][is1][is4]*GOS[imu5][ic3][ic2][is3][is2]-GTW[imu5][ic1][ic4][is1][is4]*GOS[imu][ic3][ic2][is3][is2]);
        }
        // SP -PS
                OpmPV[2][ic1][ic2][ic3][ic4][is1][is2][is3][is4]=OpmPV[2][ic1][ic2][ic3][ic4][is1][is2][is3][is4]+
                1./2.*(GOS[0][ic1][ic2][is1][is2]*GTW[5][ic3][ic4][is3][is4]-GOS[5][ic1][ic2][is1][is2]*GTW[0][ic3][ic4][is3][is4])
                -pm(ipm)*1./2.*(GTW[0][ic1][ic4][is1][is4]*GOS[5][ic3][ic2][is3][is2]-GTW[5][ic1][ic4][is1][is4]*GOS[0][ic3][ic2][is3][is2]);

        // SP+PS
                OpmPV[3][ic1][ic2][ic3][ic4][is1][is2][is3][is4]=OpmPV[3][ic1][ic2][ic3][ic4][is1][is2][is3][is4]+
                1./2.*(GOS[0][ic1][ic2][is1][is2]*GTW[5][ic3][ic4][is3][is4]+GOS[5][ic1][ic2][is1][is2]*GTW[0][ic3][ic4][is3][is4])
                -pm(ipm)*1./2.*(GTW[0][ic1][ic4][is1][is4]*GOS[5][ic3][ic2][is3][is2]+GTW[5][ic1][ic4][is1][is4]*GOS[0][ic3][ic2][is3][is2]);

        // TT~ 
        for (int isigma=6; isigma<12; isigma++)
        {
                int isigma5=sigma_dual(isigma);
 	        OpmPV[4][ic1][ic2][ic3][ic4][is1][is2][is3][is4]=OpmPV[4][ic1][ic2][ic3][ic4][is1][is2][is3][is4]+
                1./2.*(GOS[isigma][ic1][ic2][is1][is2]*GTW[isigma5][ic3][ic4][is3][is4])
               -pm(ipm)*1./2.*(GTW[isigma][ic1][ic4][is1][is4]*GOS[isigma5][ic3][ic2][is3][is2]);



        }


   }// end ccccssss 



}
void compute_OpPC(O5_ccccssss OpmPC_out, gamma_ccss GOS, gamma_ccss GTW)
{
        compute_OPC(OpmPC_out,GOS,GTW,0);
}

void compute_OmPC(O5_ccccssss OpmPC_out, gamma_ccss GOS, gamma_ccss GTW)
{
        compute_OPC(OpmPC_out,GOS,GTW,1);
}

void compute_OpPV(O5_ccccssss OpmPV_out, gamma_ccss GOS, gamma_ccss GTW)
{
        compute_OPV(OpmPV_out,GOS,GTW,0);
}

void compute_OmPV(O5_ccccssss OpmPV_out, gamma_ccss GOS, gamma_ccss GTW)
{
        compute_OPV(OpmPV_out,GOS,GTW,1);
}




/////////////////////////////////////////////////////////////////////////////////////
// routines computing the two relevant r combinations
// 0001 & 1110
///////////////////////////////////////////////
// Recordemos que:
// Como:
// OpmPX(0010)=OpmPX(1011)=OpmPX(1110)*=OpmPX(1000)*
// OpmPX(0001)=OpmPX(1101)=OpmPX(0111)*=OpmPX(0100)*
// (Comprobado numéricamente)
// (Parece lógico porque es cambiar la parte TW - OS pero buscar la simetria) //TODO
// Basta calcular
// OpmPX(0001) y OpmPX(1110) y al promediar quedarnos con la parte real
//////////////////////////////////////////////
// No hay suficiente memoria para correr estas rutinas que nos harian las dos r a la vez
////////////////////////////////////////////////////////////////////////////////////////

void compute_OPC_2r(r_O5_ccccssss OpmPC_out, gamma_rr_ccss G,int ipm)
{


        O5_ccccssss ****OpmPC=malloc(sizeof(O5_ccccssss***)*(2));
        for (int r1=0; r1<2; r1++) OpmPC[r1]=malloc(sizeof(O5_ccccssss**)*(2));
        for (int r1=0; r1<2; r1++) for (int r2=0; r2<2; r2++)  OpmPC[r1][r2]=malloc(sizeof(O5_ccccssss*)*(2));
        for (int r1=0; r1<2; r1++) for (int r2=0; r2<2; r2++) for (int r3=0; r3<2; r3++) OpmPC[r1][r2][r3]=malloc(sizeof(O5_ccccssss)*(2));


     // 1. Initialize
      for (int ic1=0; ic1<3 ; ic1++) for (int ic2=0; ic2<3 ; ic2++) for (int ic3=0; ic3<3 ; ic3++) for (int ic4=0; ic4<3 ; ic4++)
      for (int is1=0; is1<4 ; is1++) for (int is2=0; is2<4 ; is2++) for (int is3=0; is3<4 ; is3++) for (int is4=0; is4<4 ; is4++)
      for (int iop=0; iop<5; iop++)
      for (int r1=0; r1<2; r1++) for (int r2=0; r2<2; r2++) for (int r3=0; r3<2; r3++) for (int r4=0; r4<2; r4++)
      {
                OpmPC[r1][r2][r3][r4][iop][ic1][ic2][ic3][ic4][is1][is2][is3][is4]=0;
      }

        // 2. Compute for each r combination with one r different

      for (int ircomb=0; ircomb<2 ; ircomb++)
      {
         // Compute only r-combination of r1=r2=r3=-r4 type
         int *rvector=malloc(sizeof(int)*(4)); whichr_sort(ircomb,rvector);
         int r1,r2,r3,r4; r1=rvector[0]; r2=rvector[1]; r3=rvector[2]; r4=rvector[3];

        for (int ic1=0; ic1<3 ; ic1++) for (int ic2=0; ic2<3 ; ic2++) for (int ic3=0; ic3<3 ; ic3++) for (int ic4=0; ic4<3 ; ic4++)
        for (int is1=0; is1<4 ; is1++) for (int is2=0; is2<4 ; is2++) for (int is3=0; is3<4 ; is3++) for (int is4=0; is4<4 ; is4++)
        {


         // VV+AA

         for (int imu=1; imu<5 ; imu++)
         {
                int imu5=imu+11;

                OpmPC[r1][r2][r3][r4][0][ic1][ic2][ic3][ic4][is1][is2][is3][is4]=OpmPC[r1][r2][r3][r4][0][ic1][ic2][ic3][ic4][is1][is2][is3][is4]+
                1./2.*(G[imu][r1][r2][ic1][ic2][is1][is2]*G[imu][r3][r4][ic3][ic4][is3][is4]+G[imu5][r1][r2][ic1][ic2][is1][is2]*G[imu5][r3][r4][ic3][ic4][is3][is4])
                -pm(ipm)*1./2.*(G[imu][r1][r4][ic1][ic4][is1][is4]*G[imu][r3][r2][ic3][ic2][is3][is2]+G[imu5][r1][r4][ic1][ic4][is1][is4]*G[imu5][r3][r2][ic3][ic2][is3][is2]);
        }
         // VV-AA

         for (int imu=1; imu<5 ; imu++)
         {
                int imu5=imu+11;
                OpmPC[r1][r2][r3][r4][1][ic1][ic2][ic3][ic4][is1][is2][is3][is4]=OpmPC[r1][r2][r3][r4][1][ic1][ic2][ic3][ic4][is1][is2][is3][is4]+
                1./2.*(G[imu][r1][r2][ic1][ic2][is1][is2]*G[imu][r3][r4][ic3][ic4][is3][is4]-G[imu5][r1][r2][ic1][ic2][is1][is2]*G[imu5][r3][r4][ic3][ic4][is3][is4])
                -pm(ipm)*1./2.*(G[imu][r1][r4][ic1][ic4][is1][is4]*G[imu][r3][r2][ic3][ic2][is3][is2]-G[imu5][r1][r4][ic1][ic4][is1][is4]*G[imu5][r3][r2][ic3][ic2][is3][is2]);
        }

        // SS - PP
                OpmPC[r1][r2][r3][r4][2][ic1][ic2][ic3][ic4][is1][is2][is3][is4]=OpmPC[r1][r2][r3][r4][2][ic1][ic2][ic3][ic4][is1][is2][is3][is4]+
                1./2.*(G[0][r1][r2][ic1][ic2][is1][is2]*G[0][r3][r4][ic3][ic4][is3][is4]-G[5][r1][r2][ic1][ic2][is1][is2]*G[5][r3][r4][ic3][ic4][is3][is4])
                -pm(ipm)*1./2.*(G[0][r1][r4][ic1][ic4][is1][is4]*G[0][r3][r2][ic3][ic2][is3][is2]-G[5][r1][r4][ic1][ic4][is1][is4]*G[5][r3][r2][ic3][ic2][is3][is2]);


         // SS + PP
                OpmPC[r1][r2][r3][r4][3][ic1][ic2][ic3][ic4][is1][is2][is3][is4]=OpmPC[r1][r2][r3][r4][3][ic1][ic2][ic3][ic4][is1][is2][is3][is4]+
                1./2.*(G[0][r1][r2][ic1][ic2][is1][is2]*G[0][r3][r4][ic3][ic4][is3][is4]+G[5][r1][r2][ic1][ic2][is1][is2]*G[5][r3][r4][ic3][ic4][is3][is4])
                -pm(ipm)*1./2.*(G[0][r1][r4][ic1][ic4][is1][is4]*G[0][r3][r2][ic3][ic2][is3][is2]+G[5][r1][r4][ic1][ic4][is1][is4]*G[5][r3][r2][ic3][ic2][is3][is2]);

        // TT
        for (int isigma=6; isigma<12; isigma++)
        {
                OpmPC[r1][r2][r3][r4][4][ic1][ic2][ic3][ic4][is1][is2][is3][is4]=OpmPC[r1][r2][r3][r4][4][ic1][ic2][ic3][ic4][is1][is2][is3][is4]+
                1./2.*(G[isigma][r1][r2][ic1][ic2][is1][is2]*G[isigma][r3][r4][ic3][ic4][is3][is4])
                -pm(ipm)*1./2.*(G[isigma][r1][r4][ic1][ic4][is1][is4]*G[isigma][r3][r2][ic3][ic2][is3][is2]);

        }

        for (int iop=0; iop<5; iop++)
          OpmPC_out[ircomb][iop][ic1][ic2][ic3][ic4][is1][is2][is3][is4]=OpmPC[r1][r2][r3][r4][iop][ic1][ic2][ic3][ic4][is1][is2][is3][is4];

        }// end ccccssss 

   } // end r-combination

}
void compute_OPV_2r(r_O5_ccccssss OpmPV_out, gamma_rr_ccss G, int ipm)
{
        O5_ccccssss ****OpmPV=malloc(sizeof(O5_ccccssss***)*(2));
        for (int r1=0; r1<2; r1++) OpmPV[r1]=malloc(sizeof(O5_ccccssss**)*(2));
        for (int r1=0; r1<2; r1++) for (int r2=0; r2<2; r2++)  OpmPV[r1][r2]=malloc(sizeof(O5_ccccssss*)*(2));
        for (int r1=0; r1<2; r1++) for (int r2=0; r2<2; r2++) for (int r3=0; r3<2; r3++) OpmPV[r1][r2][r3]=malloc(sizeof(O5_ccccssss)*(2));


     // 1. Initialize
      for (int ic1=0; ic1<3 ; ic1++) for (int ic2=0; ic2<3 ; ic2++) for (int ic3=0; ic3<3 ; ic3++) for (int ic4=0; ic4<3 ; ic4++)
      for (int is1=0; is1<4 ; is1++) for (int is2=0; is2<4 ; is2++) for (int is3=0; is3<4 ; is3++) for (int is4=0; is4<4 ; is4++)
      for (int iop=0; iop<5; iop++)
      for (int r1=0; r1<2; r1++) for (int r2=0; r2<2; r2++) for (int r3=0; r3<2; r3++) for (int r4=0; r4<2; r4++)

      {
             OpmPV[r1][r2][r3][r4][iop][ic1][ic2][ic3][ic4][is1][is2][is3][is4]=0;
      }


     for (int ircomb=0; ircomb<2 ; ircomb++)
     {
         int *rvector=malloc(sizeof(int)*(4)); whichr_sort(ircomb,rvector);
         int r1,r2,r3,r4; r1=rvector[0]; r2=rvector[1]; r3=rvector[2]; r4=rvector[3];


        for (int ic1=0; ic1<3 ; ic1++) for (int ic2=0; ic2<3 ; ic2++) for (int ic3=0; ic3<3 ; ic3++) for (int ic4=0; ic4<3 ; ic4++)
        for (int is1=0; is1<4 ; is1++) for (int is2=0; is2<4 ; is2++) for (int is3=0; is3<4 ; is3++) for (int is4=0; is4<4 ; is4++)

        {
        // VA+AV
         for (int imu=1; imu<5 ; imu++)
         {
                int imu5=imu+11;
                OpmPV[r1][r2][r3][r4][0][ic1][ic2][ic3][ic4][is1][is2][is3][is4]=OpmPV[r1][r2][r3][r4][0][ic1][ic2][ic3][ic4][is1][is2][is3][is4]+
                1./2.*(G[imu][r1][r2][ic1][ic2][is1][is2]*G[imu5][r3][r4][ic3][ic4][is3][is4]+G[imu5][r1][r2][ic1][ic2][is1][is2]*G[imu][r3][r4][ic3][ic4][is3][is4])
                -pm(ipm)*1./2.*(G[imu][r1][r4][ic1][ic4][is1][is4]*G[imu5][r3][r2][ic3][ic2][is3][is2]+G[imu5][r1][r4][ic1][ic4][is1][is4]*G[imu][r3][r2][ic3][ic2][is3][is2]);
        }

        // VA-AV
         for (int imu=1; imu<5 ; imu++)
         {
                int imu5=imu+11;
                OpmPV[r1][r2][r3][r4][1][ic1][ic2][ic3][ic4][is1][is2][is3][is4]=OpmPV[r1][r2][r3][r4][1][ic1][ic2][ic3][ic4][is1][is2][is3][is4]+
                1./2.*(G[imu][r1][r2][ic1][ic2][is1][is2]*G[imu5][r3][r4][ic3][ic4][is3][is4]-G[imu5][r1][r2][ic1][ic2][is1][is2]*G[imu][r3][r4][ic3][ic4][is3][is4])
                -pm(ipm)*1./2.*(G[imu][r1][r4][ic1][ic4][is1][is4]*G[imu5][r3][r2][ic3][ic2][is3][is2]-G[imu5][r1][r4][ic1][ic4][is1][is4]*G[imu][r3][r2][ic3][ic2][is3][is2]);
        }
        // SP -PS
                OpmPV[r1][r2][r3][r4][2][ic1][ic2][ic3][ic4][is1][is2][is3][is4]=OpmPV[r1][r2][r3][r4][2][ic1][ic2][ic3][ic4][is1][is2][is3][is4]+
                1./2.*(G[0][r1][r2][ic1][ic2][is1][is2]*G[5][r3][r4][ic3][ic4][is3][is4]-G[5][r1][r2][ic1][ic2][is1][is2]*G[0][r3][r4][ic3][ic4][is3][is4])
                -pm(ipm)*1./2.*(G[0][r1][r4][ic1][ic4][is1][is4]*G[5][r3][r2][ic3][ic2][is3][is2]-G[5][r1][r4][ic1][ic4][is1][is4]*G[0][r3][r2][ic3][ic2][is3][is2]);

        // SP+PS
                OpmPV[r1][r2][r3][r4][3][ic1][ic2][ic3][ic4][is1][is2][is3][is4]=OpmPV[r1][r2][r3][r4][3][ic1][ic2][ic3][ic4][is1][is2][is3][is4]+
                1./2.*(G[0][r1][r2][ic1][ic2][is1][is2]*G[5][r3][r4][ic3][ic4][is3][is4]+G[5][r1][r2][ic1][ic2][is1][is2]*G[0][r3][r4][ic3][ic4][is3][is4])
                -pm(ipm)*1./2.*(G[0][r1][r4][ic1][ic4][is1][is4]*G[5][r3][r2][ic3][ic2][is3][is2]+G[5][r1][r4][ic1][ic4][is1][is4]*G[0][r3][r2][ic3][ic2][is3][is2]);

        // TT~ 
        for (int isigma=6; isigma<12; isigma++)
        {
                int isigma5=sigma_dual(isigma);
                OpmPV[r1][r2][r3][r4][4][ic1][ic2][ic3][ic4][is1][is2][is3][is4]=OpmPV[r1][r2][r3][r4][4][ic1][ic2][ic3][ic4][is1][is2][is3][is4]+
                1./2.*(G[isigma][r1][r2][ic1][ic2][is1][is2]*G[isigma5][r3][r4][ic3][ic4][is3][is4])
               -pm(ipm)*1./2.*(G[isigma][r1][r4][ic1][ic4][is1][is4]*G[isigma5][r3][r2][ic3][ic2][is3][is2]);

        }

        for (int iop=0; iop<5; iop++)
          OpmPV_out[ircomb][iop][ic1][ic2][ic3][ic4][is1][is2][is3][is4]=OpmPV[r1][r2][r3][r4][iop][ic1][ic2][ic3][ic4][is1][is2][is3][is4];

   }// end ccccssss 


  } // end for over rcomb

}

void compute_OpPC_2r(r_O5_ccccssss OpmPC_out, gamma_rr_ccss G)
{
        compute_OPC_2r(OpmPC_out,G,0);
}

void compute_OmPC_2r(r_O5_ccccssss OpmPC_out, gamma_rr_ccss G)
{
        compute_OPC_2r(OpmPC_out,G,1);
}

void compute_OpPV_2r(r_O5_ccccssss OpmPV_out, gamma_rr_ccss G)
{
        compute_OPV_2r(OpmPV_out,G,0);
}

void compute_OmPV_2r(r_O5_ccccssss OpmPV_out, gamma_rr_ccss G)
{
        compute_OPV_2r(OpmPV_out,G,1);
}

